#include "stdafx.h"

#include "config.h"

namespace DeepZoomImage 
{
	GArray * double_garray_from_str (const char *s, gint max)
	{
		char *next = (char *)s;
		GArray *values = g_array_sized_new (false, true, sizeof (double), max > 0 ? max : 16);
		double coord = 0.0;
		guint end = max > 0 ? max : G_MAXINT;

		while (next && values->len < end) {
			while (g_ascii_isspace (*next) || *next == ',')
				next = g_utf8_next_char (next);

			if (next) {
				errno = 0;
				char *prev = next;
				coord = g_ascii_strtod (prev, &next);
				if (errno != 0 || next == prev)
					goto error;

				g_array_append_val (values, coord);
			}
		}

error:
		while (values->len < (guint) max) {
			coord = 0.0;
			g_array_append_val (values, coord);
		}

		return values;
	}

	float strtof (const char *string, char **end_ptr)
	{
		return strtod( string, end_ptr );
	}

	static void msort (void *array, void *buf, size_t low, size_t high, size_t size, int copy_mode, GCompareFunc compare)
	{
		char *a1, *al, *am, *ah, *ls, *hs, *lo, *hi, *b;
		size_t copied = 0;
		size_t mid;

		mid = MID (low, high);

		if (mid + 1 < high)
			msort (array, buf, mid + 1, high, size, copy_mode, compare);

		if (mid > low)
			msort (array, buf, low, mid, size, copy_mode, compare);

		ah = ((char *) array) + ((high + 1) * size);
		am = ((char *) array) + ((mid + 1) * size);
		a1 = al = ((char *) array) + (low * size);

		b = (char *) buf;
		lo = al;
		hi = am;

		do {
			ls = lo;
			hs = hi;

			if (lo > al || hi > am) {
				/* our last loop already compared lo & hi and found lo <= hi */
				lo += size;
			}

			while (lo < am && compare (lo, hi) <= 0)
				lo += size;

			if (lo < am) {
				if (copied == 0) {
					/* avoid copying the leading items */
					a1 = lo;
					ls = lo;
				}

				/* our last compare tells us hi < lo */
				hi += size;

				while (hi < ah && compare (hi, lo) < 0)
					hi += size;

				if (lo > ls) {
					MEMCOPY (b, ls, lo - ls);
					copied += (lo - ls);
					b += (lo - ls);
				}

				MEMCOPY (b, hs, hi - hs);
				copied += (hi - hs);
				b += (hi - hs);
			} else if (copied) {
				MEMCOPY (b, ls, lo - ls);
				copied += (lo - ls);
				b += (lo - ls);

				/* copy everything we needed to re-order back into array */
				MEMCOPY (a1, buf, copied);
				return;
			} else {
				/* everything already in order */
				return;
			}
		} while (hi < ah);

		if (lo < am) {
			MEMCOPY (b, lo, am - lo);
			copied += (am - lo);
		}

		MEMCOPY (a1, buf, copied);
	}

	void MergeSort (void *base, size_t nmemb, size_t size, GCompareFunc compare)
	{
		int copy_mode;
		void *tmp;

		if (nmemb < 2)
			return;

		if (!(tmp = malloc (nmemb * size))) {
			/* fall back to using qsort() and hope & pray that the results are "stable" */
			qsort (base, nmemb, size, compare);
			return;
		}

		if ((((char *) base) - ((char *) 0)) % sizeof (long) == 0 && (size % sizeof (long)) == 0)
			copy_mode = 1;
		else if ((((char *) base) - ((char *) 0)) % sizeof (int) == 0 && (size % sizeof (int)) == 0)
			copy_mode = 2;
		else
			copy_mode = 0;

		msort (base, tmp, 0, nmemb - 1, size, copy_mode, compare);

		free (tmp);
	}

	#define UINTTRYPARSE(bits, max)                                         \
		bool                                                                    \
		UInt##bits##TryParse (const char *str, guint##bits *retval, int *err)   \
		{                                                                       \
		register const char *inptr = str;                               \
		guint##bits val = 0;                                            \
		guint digit;                                                    \
		\
		while (*inptr == ' ')                                           \
		inptr++;                                                \
		\
		str = inptr;                                                    \
		while (*inptr >= '0' && *inptr <= '9') {                        \
		digit = (*inptr - '0');                                 \
		if (val > (max / 10)) {                                 \
		*err = EOVERFLOW;                               \
		return false;                                   \
		} else if (val == (max / 10) && digit > (max % 10)) {   \
		*err = EOVERFLOW;                               \
		return false;                                   \
		} else {                                                \
		val = (val * 10) + digit;                       \
		}                                                       \
		\
		inptr++;                                                \
		}                                                               \
		\
		while (*inptr == ' ')                                           \
		inptr++;                                                \
		\
		if (inptr == str || *inptr) {                                   \
		*err = EINVAL;                                          \
		return false;                                           \
		}                                                               \
		\
		*retval = val;                                                  \
		\
		return true;                                                    \
		}

	//UINTTRYPARSE(8, 255)
	//UINTTRYPARSE(16, 65535)
	UINTTRYPARSE(32, 4294967295UL)
	//UINTTRYPARSE(64, 18446744073709551615ULL)

	#define INTTRYPARSE(bits, max)                                          \
		bool                                                                    \
		Int##bits##TryParse (const char *str, gint##bits *retval, int *err)     \
		{                                                                       \
		register const char *inptr = str;                               \
		gint##bits val = 0;                                             \
		int digit, sign = 1;                                            \
		\
		while (*inptr == ' ')                                           \
		inptr++;                                                \
		\
		if (*inptr == '-') {                                            \
		sign = -1;                                              \
		inptr++;                                                \
		}                                                               \
		\
		str = inptr;                                                    \
		while (*inptr >= '0' && *inptr <= '9') {                        \
		digit = (*inptr - '0');                                 \
		if (val > (max / 10)) {                                 \
		*err = EOVERFLOW;                               \
		return false;                                   \
		} else if (val == (max / 10)) {                         \
		if (digit > (max % 10) &&                       \
		(sign > 0 || digit > ((max % 10) + 1))) {   \
		*err = EOVERFLOW;                       \
		return false;                           \
		}                                               \
		\
		if (sign < 0)                                   \
		val = (val * sign * 10) - digit;        \
				else                                            \
				val = (val * 10) + digit;               \
				\
				inptr++;                                        \
				\
				goto done;                                      \
		} else {                                                \
		val = (val * 10) + digit;                       \
		}                                                       \
		\
		inptr++;                                                \
		}                                                               \
		\
		val *= sign;                                                    \
		\
	done:                                                                  \
		\
		while (*inptr == ' ')                                           \
		inptr++;                                                \
		\
		if (inptr == str || *inptr) {                                   \
		*err = EINVAL;                                          \
		return false;                                           \
		}                                                               \
		\
		*retval = val;                                                  \
		\
		return true;                                                    \
		}

		//INTTRYPARSE(8, 127)
		//INTTRYPARSE(16, 32767)
		INTTRYPARSE(32, 2147483647L)
		//INTTRYPARSE(64, 9223372036854775807LL)

#if 0
#define DECIMAL_DIG	15 

	float _strtof (const char *string, char **end_ptr)
	{
		const char *end;
		long double z;

		end = __atod (&z, string, FLT_MIN_EXP, FLT_MAX_EXP, 127, FLT_MANT_DIG, DECIMAL_DIG, FLT_MAX_10_EXP, -45);
		if (end_ptr != NULL)
			*end_ptr = (char *)end;
		return (float)z;
	}

#define BBITS           (LDBL_MAX_EXP + LDBL_MANT_DIG - LDBL_MAX_10_EXP \
	+ (DECIMAL_DIG * 10)/3)
#define BWORDS          (_BI_WORDS (BBITS) + 1)


#define MAX1    (_BI_WORDMAX / 10)
#define MAX2    (_BI_WORDMAX - 10 * MAX1)

#define DOT     (unsigned char)_cur_lconv.decimal_point[0]

	const char * __atod (long double *p_result, const char *string,
		int min_exp, int max_exp, int bias,
		int mant_dig, int decimal_dig,
		int max_10_exp, int min_den_10_exp)
	{
		const char *digits, *digits_end, *s, *end;
		int e, k, i, cmp, int_digits, frac_digits, ov, pt, next, neg;
		unsigned long long mant;
		_bi_word w1, w2;
		_BI_DECLARE (f, BWORDS);
		_BI_DECLARE (t, BWORDS);      /* r, temporary */
		_BI_DECLARE (x, BWORDS);      /* x and q */
		_BI_DECLARE (uv, BWORDS);     /* u or v */
		union
		{
			long double ld;
			struct
			{
				unsigned long long mant;
				unsigned short exp;
			} x;
		} z;

#define INTEGER_BIT     (1ULL << (mant_dig - 1))

		_BI_INIT (f);
		_BI_INIT (t);
		_BI_INIT (x);
		_BI_INIT (uv);

		s = string;
		while (*s != 0 && isspace ((unsigned char)*s))
			++s;
		neg = 0;
		if (*s == '-')
		{
			++s; neg = 1;
		}
		else if (*s == '+')
			++s;

		if ((s[0] == 'i' || s[0] == 'I')
			&& (s[1] == 'n' || s[1] == 'N')
			&& (s[2] == 'f' || s[2] == 'F'))
		{
			s += 3;
			if ((s[0] == 'i' || s[0] == 'I')
				&& (s[1] == 'n' || s[1] == 'N')
				&& (s[2] == 'i' || s[2] == 'I')
				&& (s[3] == 't' || s[3] == 'T')
				&& (s[4] == 'y' || s[4] == 'Y'))
				s += 5;
			*p_result = neg ? -INFINITY : INFINITY;
			return s;
		}
		else if ((s[0] == 'n' || s[0] == 'N')
			&& (s[1] == 'a' || s[1] == 'A')
			&& (s[2] == 'n' || s[2] == 'N'))
		{
			s += 3;
			end = s;
			if (end[0] == '(')
			{
				++end;
				while (*end == '_' || (*end >= '0' && *end <= '9')
					|| (*end >= 'a' && *end <= 'z')
					|| (*end >= 'A' && *end <= 'Z'))
					++end;
				if (*end == ')')
					++end;
				else
					end = s;
			}
			*p_result = NAN;
			return end;
		}

		digits = s;
		pt = 0; int_digits = frac_digits = 0;
		while ((*s >= '0' && *s <= '9') || (*s == DOT && !pt))
		{
			if (*s++ == DOT)
				pt = 1;
			else if (pt)
				++frac_digits;
			else
				++int_digits;
		}
		digits_end = s;

		if (int_digits == 0 && frac_digits == 0)
		{
			/* Subject sequence is empty or does not have the expected
			form. */

			*p_result = 0.0;
			return string;
		}

		/* Read exponent. */

		e = 0; end = s;
		if (*s == 'e' || *s == 'E')
		{
			char *exp_end, exp_neg;
			long exp;
			int saved_errno;

			++s; exp_neg = 0;
			if (*s == '+')
				++s;
			else if (*s == '-')
			{
				exp_neg = 1; ++s;
			}
			if (isdigit ((unsigned char)*s)) /* Note: strtol() skips whitespace! */
			{
				saved_errno = errno; errno = 0;
				exp = strtol (s, &exp_end, 10);
				if (errno == ERANGE)
				{
					/* If the exponent overflows, we need a *very* long
					mantissa to get a valid number.  Such a mantissa
					doesn't fit into memory. */

					end = exp_end;
					if (exp_neg)
						goto underflow;
					else
						goto overflow;
				}
				else if (errno == 0 && exp_end != s)
				{
					end = exp_end;
					e = exp_neg ? -exp : exp;
				}
				errno = saved_errno;
			}
		}

		/* Adjust the exponent so that we can treat the sequence of digits
		as integer (by ignoring the decimal point). */

		e -= frac_digits;

		/* Remove leading zeros to keep numbers small. */

		while (digits < digits_end && (*digits == '0' || *digits == DOT))
			++digits;

		/* Remove trailing zeros to keep numbers small. */

		while (digits_end > digits
			&& (digits_end[-1] == '0' || digits_end[-1] == DOT))
		{
			if (digits_end[-1] != DOT)
				++e;
			--digits_end;
		}

		/* Return 0.0 if the number is zero. */

		if (digits == digits_end)
		{
			*p_result = neg ? -0.0 : 0.0;
			return end;
		}

		/* Quick test of the exponent.  Note that the string may or may not
		contain a decimal point. */

		if (e + digits_end - digits - 2 > max_10_exp)
			goto overflow;
		if (e + digits_end - digits < min_den_10_exp)
			goto underflow;

		/* Build the integer. */

		f.n = 0;                      /* f := 0 */
		w1 = 0;                       /* Accumulates up to 9 digits */
		w2 = 1;                       /* 10 ^ number_of_digits */
		k = 0;                        /* Truncate after DECIMAL_DIG digits */
		for (s = digits; s < digits_end; ++s)
			if (*s != DOT)
			{
				if (k == decimal_dig)
					++e;
				else
				{
					/* We have to check both W1 and W2 for overflow. */

					if (w2 >= MAX1 || (w1 >= MAX1 && (w1 > MAX1 || *s - '0' > MAX2)))
					{
						/* Overflow of the accumulator.  Append the digits
						from W1 to F. */

						ov = _bi_mul_bw (&f, BWORDS, &f, w2);
						assert (!ov);
						if (w1 != 0)
						{
							t.n = 1; t.v[0] = w1;
							ov = _bi_add_bb (&f, BWORDS, &f, &t);
							assert (!ov);
						}
						w1 = 0; w2 = 1;
					}
					w1 = w1 * 10 + *s - '0';
					w2 = w2 * 10;
					++k;
				}
			}

			/* Check the exponent again, based on k -- now we don't have to cope
			with the decimal point and get a better result. */

			if (e + k - 1 > max_10_exp)
				goto overflow;
			if (e + k < min_den_10_exp)
				goto underflow;

			/* If there are digits accumulated in W1, append them to F. */

			if (w2 != 1)
			{
				ov = _bi_mul_bw (&f, BWORDS, &f, w2);
				assert (!ov);
				t.n = 1; t.v[0] = w1;
				ov = _bi_add_bb (&f, BWORDS, &f, &t);
				assert (!ov);
			}

			/* Now, the number to convert is F*10^E.  The following code is
			based on AlgorithmM of [1]. */

			k = bias + mant_dig - 1;

			/* We use powers of 5, not 10, to make the numbers smaller.  Adjust
			K appropriately. */

			k += e;

			if (e < 0)
			{
				ov = _bi_pow5 (&uv, BWORDS, -e, 0, NULL);
				assert (!ov);

				/* Estimate the scaling factor. */

				i = _bi_fls (&f) - _bi_fls (&uv);
				if (k + i - mant_dig < 1)
					i = mant_dig - k + 1;
				k += i - mant_dig;
				if (i < mant_dig)
				{
					ov = _bi_shl_b (&f, BWORDS, &f, mant_dig - i);
					assert (!ov);
				}
				else if (i > mant_dig)
				{
					ov = _bi_shl_b (&uv, BWORDS, &uv,  i - mant_dig);
					assert (!ov);
				}

				_bi_set_b (&t, BWORDS, &f);
				ov = _bi_div_rem_bb (&t, BWORDS, &x, BWORDS, &uv); /* x:=r/v, r:=r%v */
				assert (!ov);
				i = _bi_fls (&x);

				/* Our estimate for K is never too big, except for denormals. */

				if (i < mant_dig)
				{
					--k;
					assert (k == 0);
				}

				/* Sometimes, our estimate for K is too small. */

				else if (i > mant_dig)
				{
					++k;
					if (k > max_exp - min_exp + 1)
						goto overflow;
					ov = _bi_shl_b (&uv, BWORDS, &uv, 1);
					assert (!ov);
					_bi_set_b (&t, BWORDS, &f);

					/* We have to compute the quotient and remainder again. */

					ov = _bi_div_rem_bb (&t, BWORDS, &x, BWORDS, &uv);
					assert (!ov);

					/* Our estimate can never be off by more than one, according
					to Theorem 1. */

					assert (_bi_fls (&x) == mant_dig);
				}

				/* r < v - r  is equivalent to 2r < v. */

				ov = _bi_shl_b (&t, BWORDS, &t, 1);
				assert (!ov);
				cmp = _bi_cmp_bb (&t, &uv);
			}
			else
			{
				ov = _bi_pow5 (&t, BWORDS, e, 0, NULL);
				assert (!ov);
				ov = _bi_mul_bb (&uv, BWORDS, &f, &t);
				assert (!ov);

				/* v is a power of two, computing the scaling factor is
				trivial. */

				i = _bi_fls (&uv);
				k += i - mant_dig;
				if (i < mant_dig)
				{
					ov = _bi_shl_b (&uv, BWORDS, &uv, mant_dig - i);
					assert (!ov);
					i = 0;
				}
				else if (i > mant_dig)
					i = i - mant_dig;
				else
					i = 0;

				/* Compute x:=u/v and u:=u MOD v where v=2^i. */

				ov = _bi_div_rem_pow2 (&uv, &x, BWORDS, i);
				assert (!ov);

				assert (x.n == 2);

				/* r < v - r is equivalent to 2r < v.  Note that we have u=r and
				v=2^i here. */

				cmp = _bi_cmp_pow2 (&uv, i - 1);
			}

			if (k > max_exp - min_exp + 1)
				goto overflow;

			assert (x.n <= 2);            /* x.n can be < 2 for denormals */
			if (x.n == 0)
				mant = 0;
			else if (x.n == 1)
				mant = x.v[0];
			else
				mant = ((unsigned long long)x.v[1] << _BI_WORDSIZE) | x.v[0];

			if (cmp < 0)
				next = 0;
			else if (cmp > 0)
				next = 1;
			else
				next = x.v[0] & 1;
			if (next)
			{
				if ((mant & ~INTEGER_BIT) == INTEGER_BIT - 1)
				{
					mant = INTEGER_BIT;
					++k;
					if (k > max_exp - min_exp + 1)
						goto overflow;
				}
				else
					++mant;
			}

			if (mant == 0)
				goto underflow;

			/* Build the floating point number from the sign, exponent, and
			significant. */

			switch (mant_dig)
			{
			case 24:
				/* Clear the integer bit. */
				mant &= ~(1ULL << (FLT_MANT_DIG - 1));
				/* Insert the exponent. */
				mant |= ((unsigned long long)k << (FLT_MANT_DIG - 1));
				if (neg)
					mant |= 0x80000000ULL;
				*p_result = *(float *)&mant;
				break;

			case 53:
				/* Clear the integer bit. */
				mant &= ~(1ULL << (DBL_MANT_DIG - 1));
				/* Insert the exponent. */
				mant |= ((unsigned long long)k << (DBL_MANT_DIG - 1));
				if (neg)
					mant |= 0x8000000000000000ULL;
				*p_result = *(double *)&mant;
				break;

			case 64:
				z.x.mant = mant;
				z.x.exp = k;
				if (neg)
					z.x.exp |= 0x8000;
				*p_result = z.ld;
				break;
			}
			return end;

overflow:
			*p_result = neg ? -HUGE_VALL : HUGE_VALL;
			errno = ERANGE;
			return end;

underflow:
			*p_result = neg ? -0.0 : 0.0;
			errno = ERANGE;
			return end;
	}


	/* Theorem 1.  If u is a normalized m+n-bit number and v is a
	normalized n-bit number, then 2^(m-1) <= [u/v] < 2^(m+1).

	Proof: To get a lower bound for [u/v] consider the worst case
	of u=1000... (u=2^(m+n-1)) and v=1111... (v=2^n - 1).

	[u/v] >= [2^(m+n-1) / (2^n - 1)]
	>= [2^(m+n-1) / 2^n]
	= 2^(m-1)

	To get an upper bound for [u/v] consider the worst case of
	u=1111... (u=2^(m+n)-1) and v=1000... (v=2^(n-1)).

	[u/v] <= [(2^(m+n) - 1) / 2^(n-1)]
	= [2^(m+1) - 2^(1-n)]
	< [2^(m+1)]
	= 2^(m+1)                           qed. */

	/* References:

	[1] William D. Clinger,
	How to read floating point numbers accurately,
	ACM SIGPLAN notices Volume 25 Number 6, June 1990. */

#endif
};
